function ava = glm_Anova(hfile, factors, betas, subjects, group, covariate)
% GLM::Anova  - create an ANOVA object from a GLM
%
% FORMAT:       ava = glm.Anova(factors, betas)
%
% Input fields:
%
%       factors     1xF struct with fields:
%        .fname     factor name
%        .ftype     'between' or 'within'
%        .levels    1xL cell array with level labels
%       betas       N-dimensional array with unique beta assignments,
%                   number of dimensions = number of within factors
%       subjects    list of subjects to include (by default: all)
%       group       SxB double array with group assignment(s)
%                   (default: empty, no between factor)
%       covariate   CxS covariate data
%
% Note: currently, only the models supported by BVQX can be created:
%       - 1 within factor
%       - 2 within factors
%       - 3 within factors
%       - 1 between factor X 1 within factor

% Version:  v0.7a
% Build:    7072411
% Date:     Jul-24 2007, 11:34 AM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm') || ...
   ~isstruct(factors) || ...
    isempty(factors) || ...
    numel(factors) > 3 || ...
   ~isfield(factors, 'fname') || ...
   ~isfield(factors, 'ftype') || ...
   ~isfield(factors, 'levels') || ...
    isempty(betas) || ...
   (~isa(betas, 'double') && ...
    (~iscell(betas) || ...
     ~ischar(betas{1}) || ...
      isempty(betas{1})))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% get content
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;

% requires RFX GLM input type
if bc.ProjectTypeRFX ~= 1
    error( ...
        'BVQXfile:InvalidObject', ...
        'Only defined for RFX GLM.' ...
    );
end

% check other arguments
if nargin < 4 || ...
   ~isa(subjects, 'double') || ...
    numel(subjects) ~= length(subjects) || ...
    ndims(subjects) > 2 || ...
    any(isnan(subjects) | subjects < 1 | subjects > bc.NrOfSubjects | subjects ~= fix(subjects)) || ...
    numel(unique(subjects)) ~= numel(subjects)
    subjects = 1:bc.NrOfSubjects;
else
    subjects = subjects(:)';
end
numsubs = numel(subjects);
if nargin < 5 || ...
   ~isa(group, 'double') || ...
    size(group, 1) ~= numsubs || ...
    any(isnan(group(:)) | group(:) < 0 | group(:) >= numsubs | group(:) ~= fix(group(:)))
    group = zeros(numsubs, 0);
end
if nargin < 6 || ...
   ~isa(covariate, 'double') || ...
    size(covariate, 1) ~= numsubs || ...
    any(isinf(covariate(:)) | isnan(covariate(:)))
    covariate = zeros(numsubs, 0);
end

% get some GLM settings
prjtype = bc.ProjectType + 1;
numsprd = bc.NrOfSubjectPredictors;
numstud = numel(bc.Study);
stfiles = cell(numstud, 1);
for stc = 1:numstud
    stfiles{stc} = bc.Study(stc).NameOfAnalyzedFile;
end
mapdims = size(bc.GLMData.RFXGlobalMap);

% generate predictor structs
bfs = cell2struct(cell(0, 0, 4), ...
    {'Type', 'Name', 'NrOfLevels', 'LevelNames'}, 3);
wfs = cell2struct(cell(0, 0, 4), ...
    {'Type', 'Name', 'NrOfLevels', 'LevelNames'}, 3);

% parse factors
numfac = numel(factors);
numbf = 0;
numwf = 0;
blevs = 0;
wlevs = [1 1 1];
for fc = 1:numfac
    if ~ischar(factors(fc).ftype) || ...
       ~ischar(factors(fc).fname) || ...
        isempty(factors(fc).fname) || ...
       ~any(strcmpi(factors(fc).ftype, {'between', 'within'}))
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid factor type/name for factor %d.', ...
            fc ...
        );
    end
    flevels = factors(fc).levels;
    switch (lower(factors(fc).ftype))
        case {'between'}
            numbf = numbf + 1;
            blevs = numel(flevels);
            if size(group, 2) < numbf || ...
                numel(unique(group(:, numbf))) ~= blevs
                error( ...
                    'BVQXfile:BadArgument', ...
                    'Invalid between-subjects factor structure.' ...
                );
            end
            bfs(numbf).Type = 1;
            bfs(numbf).Name = factors(fc).fname;
            bfs(numbf).NrOfLevels = blevs;
            for lc = 1:blevs
                if ~ischar(flevels{lc}) || ...
                    isempty(flevels{lc})
                    flevels{lc} = sprintf('Level %d', lc);
                else
                    flevels{lc} = flevels{lc}(:)';
                end
            end
            bfs(numbf).LevelNames = flevels(:);
        case {'within'}
            numwf = numwf + 1;
            wlevs(numwf) = numel(flevels);
            wfs(numwf).Type = 1;
            wfs(numwf).Name = factors(fc).fname;
            wfs(numwf).NrOfLevels = wlevs(numwf);
            for lc = 1:wlevs(numwf)
                if ~ischar(flevels{lc}) || ...
                    isempty(flevels{lc})
                    flevels{lc} = sprintf('Level %d', lc);
                else
                    flevels{lc} = flevels{lc}(:)';
                end
            end
            wfs(numwf).LevelNames = flevels(:);
    end
end

% reject between factor if not exactly one within factor
if numbf > 1 || ...
   (numbf == 1 && numwf ~= 1) || ...
    numwf == 0 || ...
    numwf > 3
    error( ...
        'BVQXfile:BadArgument', ...
        'Requested model (factors combination) not yet supported.' ...
    );
end

% get subject grouping
llevel = cell(1, numbf);
for buc = 1:numbf
    blevu = unique(group(:, buc));
    blevs = cell(1, bfs(buc).NrOfLevels);
    for lc = 1:numel(blevs)
        blevs{lc} = find(group(:, buc) == blevu(lc));
    end
    llevel{buc} = blevs;
end

% determine model type
mtype = 0;
bbeta = false;
switch (numwf)
    case {1}
        if numbf == 0
            mtype = 1;
        else
            mtype = 3;
        end
        if numel(betas) ~= size(betas, 1)
            bbeta = true;
        end
    case {2}
        mtype = 2;
        if numel(size(betas)) ~= 2 || ...
            any(size(betas) ~= wlevs(1:2))
            bbeta = true;
        end
    case {3}
        mtype = 5;
        if numel(size(betas)) ~= 3 || ...
            any(size(betas) ~= wlevs)
            bbeta = true;
        end
end

% check betas
if iscell(betas)
    dbeta = zeros(size(betas));
    for buc = 1:numel(betas)
        if ~ischar(betas{buc}) || ...
            isempty(betas{buc})
            bbeta = true;
            break;
        end
        for pc = 1:(numsprd - 1)
            if ~isempty(strfind(lower(bc.Predictor(pc).Name2), ...
                    lower(betas{buc})))
                betas{buc} = pc;
                break;
            end
        end
        if ~isa(betas{buc}, 'double')
            betas{buc} = numsprd;
        end
        dbeta(buc) = betas{buc};
    end
    betas = dbeta;
end
betas = fix(real(betas));
if any(betas(:) < 1 | betas(:) >= numsprd | isnan(betas(:))) || ...
    numel(unique(betas(:))) ~= numel(betas)
    bbeta = true;
end
if bbeta
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid dimensions or content for betas.' ...
    );
end
blist = betas(:)';
nbeta = numel(blist);

% generate new AVA object
ava = BVQXfile('new:ava');
avac = bvqxfile_getcont(ava.L);
avac.ProjectType = prjtype;
avac.ModelType = mtype;
avac.NrOfBetweenFactors = numbf;
avac.NrOfWithinFactors = numwf;
avac.NrOfCovariates = 0;
if numbf > 0
    avac.NrOfBetweenFactorLevels = blevs(:)';
else
    avac.NrOfBetweenFactorLevels = zeros(1, 0);
end
avac.NrOfWithinFactorLevels = wlevs(1:numwf);
avac.NrOfSubjects = numsubs;
avac.NrOfPredictors = numsprd;
avac.BetweenFactor = bfs;
avac.WithinFactor = wfs;
avac.ContrastName = '';
avac.Covariates = {};
avac.OriginatingGLM = sc.F;
avac.NrOfStudies = numstud;
avac.TimeCourseFiles = stfiles;
avac.Resolution = bc.Resolution;

% further calculation is somewhat based on project type
switch (prjtype)
    case {2}
        avac.XStart = bc.XStart;
        avac.XEnd   = bc.XEnd;
        avac.YStart = bc.YStart;
        avac.YEnd   = bc.YEnd;
        avac.ZStart = bc.ZStart;
        avac.ZEnd   = bc.ZEnd;
        Sbet = struct('type', '()', 'subs', ...
            {{1:mapdims(1), 1:mapdims(2), 1, 1:nbeta, 1:numsubs}});
        Smap = struct('type', '()', 'subs', {{':', ':', 1, blist}});
        Sdim = 3;
        Smin = 1;
        Sstp = 1;
        Smax = mapdims(3);
        tms = [1, 1, 1, nbeta, numsubs];
        rms = tms;
        rms(4) = 1;
    case {3}
        avac.NrOfVertices = bc.NrOfVertices;
        Sbet = struct('type', '()', 'subs', ...
            {{1:4096, 1:nbeta, 1:numsubs}});
        Smap = struct('type', '()', 'subs', {{1:4096, blist}});
        Ssmap = struct('type', '()', 'subs', {{1:4096, 1}});
        Sdim = 1;
        Smin = 1;
        Sstp = 4096;
        Smax = mapdims(1);
        tms = [1, nbeta, numsubs];
        rms = tms;
        rms(2) = 1;
    otherwise
        error( ...
            'BVQXfile:BadObject', ...
            'Only VTC and MTC GLMs supported for now.' ...
        );
end

% create necessary output maps and design matrix
sval = single(0);
desmtx = zeros(nbeta * numsubs, numsubs + nbeta + size(group, 2) + size(covariate, 2));
switch (mtype)

    % one within factor
    case {1}
        map_ms_a = repmat(sval, mapdims);
        map_ms_erra = repmat(sval, mapdims);
        for suc = 1:numsubs
            desmtx(((suc-1) * nbeta + 1):(suc * nbeta), suc) = 1;
        end
        for btc = 1:nbeta
            desmtx(btc:nbeta:end, suc + btc) = 1;
        end
        for grc = 1:size(group, 2)
            % desmtx();
        end
        for coc = 1:size(covariate, 2)
            for suc = 1:numsubs
                desmtx(((suc-1) * nbeta + 1):(suc * nbeta), end + 1 - coc) = ...
                    covariate(suc, coc);
            end
        end
        
    % two within factors
    case {2}
        map_ms_a = repmat(sval, mapdims);
        map_ms_b = repmat(sval, mapdims);
        map_ms_axb = repmat(sval, mapdims);
        map_ms_axs = repmat(sval, mapdims);
        map_ms_bxs = repmat(sval, mapdims);
        map_ms_axbxs = repmat(sval, mapdims);

    % three within factors
    case {5}
        map_ms_a = repmat(sval, mapdims);
        map_ms_b = repmat(sval, mapdims);
        map_ms_c = repmat(sval, mapdims);
        map_ms_axb = repmat(sval, mapdims);
        map_ms_axc = repmat(sval, mapdims);
        map_ms_bxc = repmat(sval, mapdims);
        map_ms_axbxc = repmat(sval, mapdims);
        map_ms_axs = repmat(sval, mapdims);
        map_ms_bxs = repmat(sval, mapdims);
        map_ms_cxs = repmat(sval, mapdims);
        map_ms_axbxs = repmat(sval, mapdims);
        map_ms_axcxs = repmat(sval, mapdims);
        map_ms_bxcxs = repmat(sval, mapdims);
        map_ms_axbxcxs = repmat(sval, mapdims);

    % one between and one within factor
    case {3}
        map_ms_a = repmat(sval, mapdims);
        map_ms_errb = repmat(sval, mapdims);
        map_ms_b = repmat(sval, mapdims);
        map_ms_axb = repmat(sval, mapdims);
        map_ms_erraxb = repmat(sval, mapdims);

end

if mtype ~= 4
    cmeans = repmat(sval, [mapdims, nbeta]);
end

% get handle to output data
glmd = bc.GLMData;

% iterate over either z dim or vertices
for Scnt = Smin:Sstp:Smax
    Sdto = min(Scnt + Sstp - 1, Smax);
    Smap.subs{Sdim} = Scnt:Sdto;
    Ssmap.subs{Sdim} = Scnt:Sdto;
    Sbet.subs{Sdim} = 1:numel(Smap.subs{Sdim});
    bmap = 0;
    bmap = subsasgn(bmap, Sbet, bmap(1));
    
    % get betas over subjects
    for suc = 1:numsubs
        Sbet.subs{end} = suc;
        bmap = subsasgn(bmap, Sbet, subsref(glmd.Subject(suc).BetaMaps, Smap));
    end
    sbmap = size(bmap);
    bmap = reshape(bmap, sbmap(1:Sdim), nbeta * numsubs)';
    
    % compute regression
    [reg{1:4}] = calcbetas(desmtx, bmap);
    
    switch (mtype)

        % one within factor
        case {1}
            m_t = mean(bmap, 1);
            % ms_t = sum((bmap - repmat(m_t, nbeta * numsubs, 1)) .^ 2)';
            m_a = (numsubs / (nbeta - 1)) * ...
                sum((reg{1}(:, numsubs+1:numsubs+nbeta) - ...
                repmat(mean(reg{1}(:, numsubs+1:numsubs+nbeta), 2), ...
                [1, nbeta])) .^ 2, 2);
            map_ms_a = subsasgn(map_ms_a, Ssmap, m_a);

        % two within factors
        case {2}
            error('not yet supported');

        % three within factors
        case {5}
            error('not yet supported');

        % one between and one within factor
        case {3}
            error('not yet supported');

    end
end

% put into AVA
switch (mtype)

    % one within factor
    case {1}
        avac.Maps.MS_A = map_ms_a;
        avac.Maps.MS_ErrorA = map_ms_erra;
        
    % two within factors
    case {2}

    % three within factors
    case {5}

    % one between and one within factor
    case {3}

end
clear map_ms_*;

if mtype ~= 4
    avac.Maps.CellMeans = cmeans;
end
clear cmeans;

% put back in storage
bvqxfile_setcont(ava.L, avac);
